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ABSTRACT 

Context. Turbulence and angular momentum transport in accretion disks remains a topic of debate. With the realization 
that dead zones are robust features of protoplanetary disks, the search for hydrodynamical sources of turbulence contin- 
ues. A possible source is the baroclinic instability (BI), which has been shown to exist in unmagnetized non-barotropic 
disks. 

Aims. We aim to verify the existence of the baroclinic instability in 3D magnetized disks, as well as its interplay with other 
instabilities, namely the magneto-rotational instability (MRI) and the magneto-elliptical instability. 

Methods. We performed local simulations of non-isothermal accretion disks with the PENCIL CODE. The entropy gradient 
that generates the baroclinic instability is linearized and included in the momentum and energy equations in the shearing 
box approximation. The model is compressible, so excitation of spiral density waves is allowed and angular momentum 
transport can be measured. 

Results. We find that the vortices generated and sustained by the baroclinic instability in the purely hydrodynamical 
regime do not survive when magnetic fields are included. The MRI by far supersedes the BI in growth rate and strength 
at saturation. The resulting turbulence is virtually identical to an MRI-only scenario. We measured the intrinsic vorticity 
profile of the vortex, finding little radial variation in the vortex core. Nevertheless, the core is disrupted by an MHD 
instability, which we identify with the magneto-elliptic instability. This instability has nearly the same range of unstable 
wavelengths as the MRI, but has higher growth rates. In fact, we identify the MRI as a limiting case of the magneto-elliptic 
instability, when the vortex aspect ratio tends to infinity (pure shear flow). We isolated its effect on the vortex, finding 
that a strong but unstable vertical magnetic field leads to channel flows inside the vortex, which stretch it apart. When 
the field is decreased or resistivity is used, we find that the vortex survives until the MRI develops in the box. The vortex 
is then destroyed by the strain of the surrounding turbulence. Constant azimuthal fields and zero net flux fields also lead 
to vortex destruction. Resistivity quenches both instabilities when the magnetic Reynolds number of the longest vertical 
wavelength of the box is near unity. 

Conclusions. We conclude that vortex excitation and self-sustenance by the baroclinic instability in protoplanetary disks is 
viable only in low ionization, i.e., the dead zone. Our results are thus in accordance with the layered accretion paradigm. 
A baroclinicly unstable dead zone should be characterized by the presence of large-scale vortices whose cores are ellipti- 
cally unstable, yet sustained by the baroclinic feedback. Since magnetic fields destroy the vortices and the MRI outweighs 
the BI, the active layers are unmodified. 

Key words. Keywords should be given 



1. Introduction 

Turbulence is the preferred mechanism for enabling ac- 
cretion in circumstellar disks, and the magneto-rotational 
instability (MRI, Balbus & Hawley 1991) is the preferred 
route to turbulence. However, the MRI requires sufficient 
ionization since the magnetic field and the gas must be 
coupled, so it should not be expected to occur in regions 
of low ionization such as the "dead zone" (Gammie 1996, 
Turner & Drake 2009). Therefore, the search for hydrody- 
namical sources of turbulence continues, if only to pro- 
vide some residual accretion in the dead zone. A dis- 
tinct possibility is the baroclinic instability (BI; Klahr & 
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Bodenheimer 2003, Klahr 2004), the interest in which has 
been recently rekindled (Petersen et al. 2007a, Petersen et 
al. 2007b, Lesur & Papaloizou 2010). 

A baroclinic flow is one where the pressure depends on 
both density and temperature, as opposed to a barotropic 
flow where the pressure only depends on density. In 
such a flow, the non-axyssimmetric misalignment be- 
tween surfaces of constant density p (isopycnals) and sur- 
faces of constant pressure p (isobars) generates vorticity. 
Mathematically, this translates into a non-zero baroclinic 
vector, V x (— p _1 Vp) = p _2 Vp x Vp. Baroclinicity has 
long been known in atmospheric dynamics to be responsi- 
ble for turbulent patterns on planets and for weather pat- 
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terns of Rossby waves (planetary waves), cyclones, and 
anticyclones on Earth. 

The difference between the baroclinic instability of 
weather patterns on planetary atmospheres and the baro- 
clinic instability in accretion disks is that the former is 
linear, whereas the latter is nonlinear (Klahr 2004, Lesur 
& Papaloizou 2010). This is because in accretion disks, 
the disturbances have to overcome the strong Keplerian 
shear that causes perturbations to be heavily dominated 
by restoring forces in all Reynolds numbers. 

The nature of the instability was clarified in the work 
of Petersen et al. (2007ab), who highlighted the importance 
of finite thermal inertia. When the thermal time is compa- 
rable to the eddy turnover time, the vortex is able to estab- 
lish an entropy gradient around itself that compensates the 
large-scale entropy gradient that created it. This entropy 
gradient back reacts on the eddy, generating more vorticity 
via buoyancy. This in turn reinforces the gradient. A pos- 
itive feedback has been established, and the eddy grows. 
This, in a nutshell, is the baroclinic instability: the sole re- 
sult of an eddy trying to counter the background entropy 
gradient that established it, and reinforcing itself by doing 
so. 

The 3D properties of the instability have been stud- 
ied by Lesur & Papaloizou (2010). They find that the vor- 
tices produced are not destroyed when baroclinicity is 
present, although they are unstable to the elliptical insta- 
bility (Kerswell 2002, Lesur & Papaloizou 2009). The satu- 
rated state of the instability is dominated by the presence 
of large-scale 3D, self-sustained, vortices with weakly tur- 
bulent cores. The study of Petersen et al. (2007ab) and most 
of that of Lesur & Papaloizou (2010) was done with spec- 
tral codes, which filter sound waves. Vortices, however, 
have the interesting property of radiating inertial-acoustic 
waves, which are known to transport angular momen- 
tum (Heinemann & Papaloizou 2009). Lesur & Papaloizou 
(2010) performed a compressible, yet 2D, simulation, with 
a resulting Shakura-Sunyaev-like oc value of 10 -3 . 

These results are intriguing, and a major question to 
ask is what their significance is in the context of the layered 
accretion paradigm. Vortices have been described in the 
literature as devoid of radial shear (Klahr & Bodenheimer 
2006), so in principle they could form and survive in the 
midst of MRI turbulence, as the simulations of Fromang 
& Nelson (2005) suggest. Moreover, if the baroclinic insta- 
bility is able to produce and sustain vortices when mag- 
netization is present, synergy with the MRI is an interest- 
ing possibility, potentially leading to higher accretion rates 
than hitherto achieved in previous works. On the other 
hand, elliptical streamlines can be heavily destabilized 
by magnetic fields (Lebovitz & Zweibel 2004, Mizerski & 
Bajer 2009). This magneto-elliptic instability may either be 
stabilized by baroclinicity, as the elliptic instability was 
shown to be (Lesur & Papaloizou 2010), or completely 
break the vortices apart thus rendering the baroclinic in- 
stability meaningless in the presence of magnetization. We 
address these open questions in this work. 

The paper is structured as follows. In Sect. 2 we in- 
troduce the model equations of the compressible shearing 
box, modified to include the contribution from the large- 
scale background entropy gradient. In view of the contro- 
versy aroused by the baroclinic instability in the literature, 
it was found prudent to establish the reliability of the nu- 
merics, as well as to provide an independent confirmation 



of the 2D results. This is done in Sect. 3. The 3D results are 
presented in Sect. 4. Table 1 and Table 2 contain summaries 
of the simulations performed for this study, referring to 
the sections and figures they are described. Conclusions 
are given in Sect. 5. 

2. The model 

We model a local patch of the disk following the shear- 
ing box approximation. The reader is referred to Regev & 
Umurhan (2008) for an extensive discussion of the limita- 
tions of the approximation. To include the baroclinic term, 
we consider a large-scale radial pressure gradient follow- 
ing a power law of index £ 

p(r) = p (r/R )~ f / 

where r is the cylindrical radius and Rq is a reference 
radius. The overbar indicates that this quantity is time- 
independent. The total pressure is ptot=P + V' wnere V ^ s 
the local fluctuation. The linearization of this gradient is 
done in the same way as the large-scale Keplerian shear is 
linearized in the shearing box. This leads to extra terms 
in the equations involving the radial pressure gradient. 
We quote the modified shearing box equations below. A 
derivation of the extra terms is presented in Appendix A. 

^ + [u-V)p = -pV-u+f D {p) (1) 

— + (w-V)M=--Vp+^ 2O (zxu) 

+ 3 -Q,u x y + §° (- - 1) x + f v (u,p) (2) 
J) A 3 

— - -n A y x + u x B-fioVJ + fM) ( 3 ) 

g + ( „.v )s ^{v, K vr>-^<I^> 

In the equations above, u is the velocity, A the mag- 
netic vector potential, B= V x A the magnetic field, /= V x 
B/}io the current (where jiQ is the magnetic permittivity), 
rj is the resistivity, T the temperature, s the entropy, and K 
the radiative conductivity. The operator 

represents the Keplerian derivative of a fluid parcel. It 

is the only place where the Keplerian flow u^p appears 
explicitly. The advection is made Galilean-invariant by 
means of the SAFI algorithm (Johansen et al. 2009), which 
speeds up performance. The simulations were done with 
the Pencil Code 1 . 

We work with entropy as the main thermodynamical 
quantity. This is a natural choice when dealing with baro- 
clinicity. Considering the polytropic equation of state, 

p = fc^ n+1 ) /n , (6) 
1 See http://www.nordita.org/software/pencil-code/ 
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Fig. 1. Snapshots of the fiducial 2D run with £ = 2, r = In, H — 0.1, and resolution 256 2 . A vortex is formed and establishes a local 
entropy gradient that counteracts the global entropy gradient that caused it in first place. Moderate cooling times keep the surfaces 
of constant density and constant pressure misaligned, leading to more vortex growth. In the positive feedback that ensues, giant 
anti-cyclonic vortices grow to the sonic scale. The initial condition was free of ens trophy. This vorticity growth was due purely to 
baroclinic effects. 
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Fig. 2. Left). Baroclinic enstrophy growth with thermal relaxation, where a gas parcel returns to the initial temperature on a timescale 
t. A stable entropy gradient can only be maintained between the extremes of too fast a relaxation (isothermal behavior) and too slow 
a relaxation (adiabatic behavior). Optimal growth occurs when r is comparable to the dynamical time. 

Right). Baroclinic enstrophy growth with thermal diffusion, where heat diffuses over a scale height on a timescale r. Optimal growth 
occurs on longer timescales when compared to the thermal relaxation case. 



and the definition of entropy, 

s = c v [ln(p/p ) - 7ln(p/p )] , (7) 

we immediately recognize s=c v ln(fc/fco) in the case 
n=l/ (7 — 1); i.e., up to a constant, entropy is the propor- 
tionality factor in the polytropic equation of state. That 
means that any spatial gradient of entropy translates into 
a departure from barotropic conditions 2 . 

The third term on the right-hand side of the entropy 
equation is an artificial thermal relaxation term, which 
drives the temperature back to the initial temperature To, 
on a pre-specified thermal timescale r c . The temperature is 

2 Actually, this is such a useful insight that some authors pre- 
fer to define entropy as S=p/p ry . The reader should then keep 
in mind that what we call entropy is actually s=c v ln(S/Srj) in 
that definition. Here we prefer to use the definition Eq. (7) as it 
comes from thermodynamics; i.e., Tds=de + pdv, where e=c v T is 
the internal energy and v=l/p. It also enables the Brunt- Vaisala 
frequency to be written in a more compact form (Eq. (8)). 



T=c 2 / [cp(y — 1)], where c s is the sound speed, y=c p /c v 
the adiabatic index, and c v and c p are the heat capacities 
at constant volume and constant pressure, respectively. 

After defining entropy, we also define the Brunt- 
VaisaTa frequency N, the frequency associated with buoy- 
ant structures 



N 2 = — —Vp-Vs (8) 

c p p 

1 dp ( 1 dp 1 dp\ 
p dr \p dr yp dr ) ' 

and we have assumed axis-symmetry (3^=0) and no ver- 
tical stratification (d z =0) between the steps. In our setup, 
there is no large-scale density gradient, so the first term 
inside the parentheses cancel. As dp/dr=-p^/r, we have, 
at r=R 
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Angular Momentum Transport 
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Fig. 3. Enstrophy and the resulting alpha-viscosity for the fidu- Fig. 4. Dependence on resolution. The low resolution run fails 
cial 2D run. The two quantities are quite well correlated, since to develop vortices, reaffirming that aliasing is not occurring in 
the angular momentum transport is the result of inertial-acoustic our models. The middle and high resolution runs saturate at the 



waves, which in turn are driven by vorticity. 



same level of enstrophy, which suggests convergence. 



(9) 



i.e., the Brunt- Vaisala frequency is always imaginary. 
However, the flow is convectively stable, since the 
epicyclic frequency squared 



r 3 dr 



(10) 



is far higher than — N 2 , so that the Solberg-Hoiland crite- 
rion is always satisfied 
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K 2 + N 2 > 0. (11) 

In Eq. (10), j=Qr 2 is the specific angular momentum per 
unit mass. 

We add explicit sixth-order hyperdif fusion /d(p), hy- 
perviscosity f v (u,p), and hyperresistivity f?j(A) to the 
mass, momentum, and induction equations as specified in 
Lyra et al (2008). Hyper heat conductivity /k(s) to the en- 
tropy equation is added as in Lyra et al. (2009) and Oishi & 
Mac Low (2009). All simulations use Cp=Ro=Qo=po=}io=l, 
7=1.4, and c s =0.1. 



Fig. 5. Dependence on Grid Reynolds number. The hyperviscosi- 
ties shown correspond to Re=0.002, 0.02, and 0.2 with respect 
to the velocity shear introduced by the Keplerian flow, calcu- 
lated on the grid scale. The initial phase of growth occurs at 
Re<l, where it is seen that the amount of growth depends on 
the Reynolds number. Upon saturation, all simulations converge 
to the same level of enstrophy. A heavily aliased solution occurs 
for v( 3 )=10 -21 , where even a simulation seeded only with noise 
develops vortices. The same does not occur for the hyperviscosi- 
ties shown, where finite amplitude perturbations were required. 
We usually use i/ 3 )=10~ 17 . 



3. 2D Results 

Although the most interesting results for the effects of the 
baroclinic generation of vorticity should be given by 3D 
models, the existence of the baroclinic instability has been 
strongly contested even in two dimensions. We therefore 
consider it important to present 2D results confirming its 
excitation. In Fig. 1 we present a fiducial 2D run where the 
evolution of the baroclinic instability is followed. It corre- 
sponds to run A in Table 1. The slope of the entropy gra- 
dient is £ = 2, which corresponds to a Richardson number 
Ri = 4N 2 /9H| ~ — 9 x 10 -3 . As initial condition, we seed 
the box with noise at small wavenumbers only, following 



Z(x,y) = CZ e- {x/2a)2 x 

E E^l^fif +;f + fc/U. (12) 
i=-k x j=o L V Lx L v /J 

The phase < (p < 1 determines the randomness. The 
subscripts underscore that the phase is the same for all 
grid points, only changing with wavenumber. The con- 
stant C sets the strength of the perturbation. As stressed by 
Lesur & Papaloizou (2010), the baroclinic instability is sub- 
critical, and therefore a finite initial amplitude is needed to 
trigger growth. We set the constant C to yield £ r ms = 0.05. 



Lyra & Klahr: Baroclinic instability in magnetized protoplanetary disks 
Table 1. Simulation suite parameters for nonmagnetic runs and results. 



5 







Parameter 


Results 






Run 


(Aorb) 


/c t h £ Resolution <<^ 2 > oc 

(/o 2 K ) 


<w 2 > 

(/<4) 


Prms 

(/po) 


<T> 
(/To) 




Fiducial 2D , Sect. 3.1, Figs. 1 and 3 


A 


1 


2 256 2 10 _iy 0.04 0.007 


0.05 


0.03 


1.01 





Thermal relaxation, Sect. 3.3, Fig. 2 



"01 2 25? KT iy 5 x 10 _i> 2 x 10" b 0.001 1~~ 

10 2 256 2 10" 17 0.01 0.001 0.01 0.01 1.01 

100 2 256 2 10" 17 6 x 10" 4 10" 5 2 x 10" 4 0.99 

oo 2 256 2 10" 17 2 x 10" 4 2 x 10" 5 0.001 0.98 



Thermal diffusion, Sect. 3.4, Fig. 2 
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0.007 


0.001 
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Resolution, Sect. 3.5, Fig. 4 

J 1 2 128x512 3 x 10" ib 7 x 10" b 2 x 10" b 1 1 

K 1 2 256x1024 10" 17 0.04 0.007 0.05 0.05 1.01 

L 1 2 512x2048 3 x 10" 18 0.04 0.008 0.06 0.05 1 - 

Reynolds number, Sect. 3.6, Fig. 5 

M 1 2 256x1024 10" ib O03 0.007 O05 O04 I - 

N 1 2 256x1024 10" 18 0.04 0.006 0.05 0.04 1.01 - 

Fiducial 3D , Sect. 4, Figs. 6, 7, 8, 13, and 16 

O I 2 256 2 xl28 lCT 17 O03 O004 O03 O02 L01 3 x 10" b 



B 
C 
D 
E 



Table 2. Simulation suite parameters for magnetic runs. 



Run t c k th g Bp \ 

MRI, Sect. 4.1, Figs. 9, 10, 11, and 16 

P 1 2 5 x 1Q- 3 

Strong field, Sect. 4.2.1, Figs. 16 and 17 

Ql 1 2 6 x IP' 1 
Q2 1 2 3.75 x IP' 2 8.8 x 10~ 4 

Weak field, Sect. 4.2.2, Figs. 16 and 18 

R 1 2 1.5 x 1Q- 3 

Resistivity, Sect. 4.2.3, Figs. 16, 19, and 20 

51 1 2 5 x 1Q- 3 1Q- 3 

52 1 2 5 x lO' 3 1.6 x IP" 4 " 

Azimuthal field, Sect. 4.3, Fig. 12 

T 1 2 3 x IP' 2 

Zero net flux field, Sect. 4.3, Fig. 12 

U 10 2 x-2 

Control runs, Sect. 4.1, Fig. 12 

CP 1 5 x lO' 3 
CR 1 1.5 x 1Q- 3 
CS 1 5 x 10~ 3 1.6 x 10" 4 



The entropy is then initialized such that p=po=const in the 
sheet. 



The rationale behind this unorthodox initial condi- 
tion is that this noise is independent of resolution. The 
usual Gaussian noise distributes power through all wave- 
lengths, so the wavelengths from k < 10 are assigned in- 
creasingly less power as the resolution increases. We stress 
that it is not vital for the instability to be seeded with 
resolution-independent noise, nor are we missing impor- 
tant physics by not exciting the small scales. 

We do not seed noise in the velocity field. The initial 
condition is strictly nonvortical. Since in 2D the stretching 
term is absent, any increase in vorticity can only be a baro- 
clinic effect. 

3.1. Baroclinic production of vorticity 
The baroclinic term, in 2D, is 

(Vp X Vp) z = d x pd y p - d x pd y p - fyoR^dyp. (13) 

The two first terms are local, whereas the third comes from 
the large-scale gradient. There would be a fourth if we con- 
sidered a large-scale density gradient as well. This third 
term generates vorticity out of any azimuthal perturba- 
tion in the density, much in the same way as the locally 
isothermal approximation does in global disks. This term 
is paramount, since it is the only source term that will gen- 
erate net enstrophy in the flow. In the beginning, this term 
dominates, generating enstrophy out of the initial density 
perturbations. The enstrophy is then amplified by the lo- 
cal baroclinic vector via the positive feedback described in 
the introduction. 
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We witness the same general phenomena as Petersen 
et al. (2007a), even though the details of implementating 
the entropy gradient are different. In global simulations, 
the vortex swings gas parcels back and forth from cold to 
hot, which causes baroclinicity and vortex growth. Here 
the initial temperature all over the box is the same, so 
the vortex does not automatically swing gas from cold 
to hot. It swings it up and down the hard-coded pres- 
sure(=temperature=entropy) gradient. The second-to-last 
term on the right-hand side of the entropy equation comes 
from the linearization of the advection term in the pres- 
ence of an entropy gradient. It embodies how the rela- 
tive entropy of a fluid parcel with respect to the back- 
ground entropy changes as it walks up or down an en- 
tropy gradient, clearly demonstrated by the dependence 
on u x . Because the movement along the vortex lines is 
embedded with a u x component of motion, this term in- 
creases or decreases the entropy within the gas parcel de- 
pending on the sign of u x . Of course, this is just the same 
physical effect in a different frame of reference. 

We show a time series of the flow in Fig. 1 as snap- 
shots of density, pressure, entropy, and vorticity. From 
these snapshots we see that the density and the pres- 
sure are very correlated. One would therefore expect that 
the amount of baroclinicity produced is tiny or vanish- 
ingly small. However, looking at the snapshots of entropy, 
we see that appearances are deceiving: the vortex gener- 
ates a strong radial entropy gradient around itself. This is 
what we described qualitatively in the introduction, and 
what Petersen et al. (2007ab) called a "sandwich pattern". 
Notice that the sign is indeed reversed with respect to the 
global gradient (higher at negative values of x). This pat- 
tern of a local entropy gradient developed by the vortex is 
a constant feature throughout the simulations. 

3.2. Angular momentum transport 

A very important question to ask is what is the strength 
of the angular momentum transport of the resulting baro- 
clinic unstable flow. We measured the kinetic alpha values 
in the simulations, 



where R x V=pSu x Suy is the Reynolds stress. The value mea- 
sured is ol « 5 x 10 -3 , indicating good transport of an- 
gular momentum. The temporal variation of alpha (Fig. 3) 
matches that of the enstrophy. 

This correlation is understood in light of the shear- 
vortex wave coupling. The angular momentum transport 
does not come from the vortex itself, but is instead caused 
by the inertial-acoustic waves that are driven by vor- 
ticity. For a detailed explanation, see Mamatsashvili & 
Chagelishvili (2007), Heinemann & Papaloizou (2009), and 
Tevzadze et al. (2010). The same production of shear waves 
and associated angular momentum transport are seen in 
the 2D compressible runs of Lesur & Papaloizou (2010). 
It should be kept in mind that the quoted angular mo- 
mentum transport may be overestimated because of us- 
ing a shearing box. As pointed out by Regev & Umurhan 
(2008), the shearing box approximation may lead to wrong 
results because it has a limited spatial scale, excessive sym- 
metries, and uses periodic boundaries. In the particular 



case of a vortex in a box, the periodic boundaries enforce 
interaction between the vortex and the strain field of its 
own images, which may lead to spurious generation of 
Reynolds stress. 

We underscore again that the initial condition was 
nonvortical. The finite-amplitude perturbations are turned 
into vortical patches by the global baroclinic term, which 
then may more owing to the local baroclinic feedback. 

Other point worth highlighting is that the instability 
has slow growth rates on the order of a hundred orbits. 
The saturated state is only weakly compressible, with the 
rms density < p 2 > — < p > 2 at a modest 0.05. 

All our simulations are compressible, so they are 
timestep limited by the presence of sound waves. The vis- 
cosity and heat diffusion are explicit, so they influence the 
Courant condition, which further limits the timestep. We 
calculated the fiducial model for 1000 orbits, but the other 
runs, which intend to explore the parameter space, were 
only calculated for 500. These are shown in the next sec- 
tions. 

3.3. Thermal time 

The fiducial simulation had a constant thermal relax- 
ation time equal to one orbit, r c =r or t, in Eq. (7), where 
r or b=27r/ 1 0. This is more representative of the very outer 
disk, ^100 AU, but at 10 AU the disk is optically thick. The 
thermal time is therefore expected to be much longer, and 
it is instructive to examine the behavior of the baroclinic 
instability in such a regime. We ran simulations with ther- 
mal times of 10 and 100 orbits, shown in Fig. 2. These cor- 
respond to runs C and D in Table 1. The extreme cases 
of an adiabatic run (run E) and a nearly isothermal one 
(r c =0.1 orbit, run B) are also presented. In agreement with 
Petersen et al. (2007a), the runs with longer thermal times 
allow for a stronger increase in enstrophy in the first or- 
bits, also seen in the adiabatic case. This is because the 
initial thermal perturbations disperse slowly without ther- 
mal relaxation, thus remaining tight (strong gradients) and 
allowing for a stronger baroclinic amplification. 

After the first eddies appear, the establishment of a 
baroclinic feedback needs a fast cooling time to lead to 
the reverted entropy gradient seen in the fiducial run. The 
most vigorous enstrophy growth in this phase is indeed 
seen to be the one with r c equal to one orbit. For a cooling 
time of 10 orbits, sustained growth of enstrophy only hap- 
pens at later times (between 20 and 50 orbits, as opposed 
to 10 orbits for t c =t ox \>), and leads to 5 times less (grid- 
averaged) enstrophy at 150 orbits. The isothermal case and 
the adiabatic cases, as expected, can never establish the 
counter entropy gradient needed for the baroclinic feed- 
back and do not experience enstrophy growth past the ini- 
tial phase. 

3.4. Diffusion 

Thermal relaxation is one of the ways of changing the in- 
ternal energy of a gas parcel. Another way is of course 
diffusion. Petersen et al. (2007a) and Lesur & Papaloizou 
(2010) report sustained baroclinic growth using thermal 
diffusion, and this is also why the 3D simulations of Klahr 
& Bodenheimer (2003) with flux-limited diffusion also ex- 
perienced baroclinic growth. In the 2D simulations by 
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Fig. 6. Evolution of vorticity (upper panels) and entropy (lower panels) due to the baroclinic instability in 3D. 



Klahr & Bodenheimer (2003), the thermal relaxation was 
numerical and a result of low resolution and a dispersive 
numerical scheme. 



We assess the effect of diffusion by setting r c =oo (i.e., 
shutting down the thermal relaxation) and adding a non- 
zero radiative conductivity K to the entropy equation (Eq. 
(7)). The thermal diffusivity k^, is related to the radiative 
conductivity by fc t h = K/ (c p p). As we choose dimension- 
less units such that Cp=po=l, then k^ « K in our simu- 
lations. The thermal diffusivity, like any diffusion coeffi- 
cient, has a dimension of L 2 T _1 , where L is length and T is 
time. We use L=H where H=c s /Q is the scale height, and 
write A: t h=H 2 /Tdiff so that the heat diffuses over one scale 
height within a time T^ff . We assess r^iff = 1, 10, and 100 or- 
bits. These are runs F, G, and H in Table 1. We see that now 
too fast a diffusion time (1 orbit) does not lead to growth, 
and vigorous growth occurs for 100 orbits. The rationale is 
the same as for the radiative case. Too fast diffusion dis- 
perses temperature gradients and weakens the baroclinic 
feedback. Slow diffusion works towards keeping the gra- 
dients tight and leads to vigorous growth. 



It is curious that the optimal diffusion time for growth 
is longer than in the thermal relaxation case. The dif- 
ference between them is that relaxation is proportional 
to the temperature, whereas diffusion is proportional to 
the Laplacian of the temperature; that is, relaxation op- 
erates equally in all spectrum, while diffusion mostly af- 
fects higher frequencies. As such, stronger diffusion (when 
compared to relaxation) should be needed for longer 
wavelengths. At present, we can offer no explanation as 
to why this is not the case, though we notice that Lesur & 
Papaloizou (2010) also see that the optimal diffusion time 
for the baroclinic feedback is substantially longer than the 
vortex turnover time. 



3.5. Resolution 

To investigate the effect of resolution, we compare runs 
using 128, 256, and 512 grid zones in the x-axis, and cells 
of unity aspect ratio (meaning four times more resolution 
in the y-axis than in the fiducial run). These are runs J, K, 
and L in Table 1. They use the default values of r c =l orbit 
for the thermal relaxation time, and £=2 for the entropy 
gradient. 

As seen in Fig. 4, the run with resolution 128x512 (low 
resolution) fails to sustain vortex growth, in contrast to 
the runs with resolution 256x1024 (middle resolution) and 
512x2048 (high resolution). That the low-resolution run 
does not lead to enstrophy growth is a salutary reassur- 
ance that aliasing is not spuriously injecting vorticity in 
the box. The high-resolution run shows a slightly higher 
initial enstrophy production (from 0-30 orbits), yet it satu- 
rates to the same level as the middle resolution run, which 
suggests convergence. 



3.6. Grid Reynolds number 

As for the grid Reynolds number, we check three runs, 
with hyper-viscosities v^=10~ 16 (run M), 10 -17 (run A), 
and 10 -18 (run N), meaning grid (hyper-)Reynolds num- 
bers of 2 x 10~ 3 , 2 x 10~ 2 , and 2 x 10 _1 with respect to 
the velocity shear introduced by the Keplerian flow, Re = 
(3/2)fi Ax 2 /v, where v=v^ (zr/Ax) 4 . 

By increasing the Reynolds number, the initial enstro- 
phy amplification is stronger. Upon saturation, the mean 
enstrophy in all simulations converge to the same value. It 
should be noticed that although the grid Reynolds num- 
ber upon saturation is greater than 1, the initial phase 
of growth occurs below this number - so growth cannot 
be due by aliasing. A heavily aliased solution is only at- 
tained when the hyperviscosity is decreased to i/( 3 )=10 -21 , 
so that the initial phase of growth occurs at very high 
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Reynolds numbers (1000). At this Reynolds number, vor- 
tex growth occurs when the simulation is seeded with 
Gaussian noise, which is a sign that the growth was nu- 
merical, given the nonlinear nature of the baroclinic in- 
stability. In contrast, none of the simulations shown in 
the figure develop vortices when only seeded with noise. 
We usually use y( 3 ">=10~ 17 which yields a good compro- 
mise between not leading to aliasing and not affecting the 
timestep too much. 



4. 3D Results 

Having examined the behavior of the baroclinic instability 
in 2D, we now turn to 3D simulations. We only study the 
unstratified case, because the stratified case needs a mod- 
ification of the evolution equations, replacing po in Eqs. 
(2) and (4) by pof(z), where f(z) is the stratification func- 
tion. We use a box of length (4xl6x2)H, with resolution 
256x256x128. Unlike in Lesur & Papaloizou (2010), our 
simulations are compressible, which limits the timestep 
and makes it impractical to follow a 3D computation for 
many hundreds of orbits. For this reason, we follow it for 
250 orbits, which was seen to be the beginning of satu- 
ration in 2D runs. The parameters of the simulation are 
shown in Table 1 as run O. 

In Fig. 6 we show snapshots of enstrophy and entropy, 
and in Fig. 7 we plot the time series of box-averaged en- 
strophy, alpha value, and rms vertical velocity. As seen 
from these figures, the 3D baroclinic instability evolves 
very similarly to its 2D counterpart. After 200 orbits the 
instability begins to saturate as vortices merge and the re- 
maining giant vortex grows to the sonic scale. The sand- 
wich pattern of entropy perturbations sustaining the vor- 
tex is also very similar. The saturated state also displays 
similar values of enstrophy (co^/fl 2 of the order of 10 -2 ) 
and angular momentum transport (oc « 5 x 10 -3 ). 

The difference is in the excitation of the elliptical insta- 
bility (Kerswell 2002, Lesur & Papaloizou 2009). As seen 
in the lower panel of Fig. 7, the growth of this instability is 
very rapid, with the rms of the vertical velocity rising by 
ten orders of magnitude in less than 10 orbits. As in Lesur 
& Papaloizou (2010), the instability leads to turbulence in 
the core of the vortex, but it is not powerful enough to 
break its coherence. This is because the elliptic destruction 
caused by the vortex stretching term is compensated with 
vorticity production by the baroclinic term. 

We follow the evolution of the vortex for 130 more 
orbits, without seeing any decay in the rms vertical ve- 
locity. In Fig. 8 we plot vertical slices of the z-vorticity 
and z-velocity, taken at the y-position of the vorticity min- 
imum, at £=250 orbits. The snapshots reveal the verti- 
cal motions at the vortex core. The motion seems turbu- 
lent, only weakly compressible, with maximum velocities 
reaching 10% of the sound speed. 

We stress again that the alpha value is around 5 x 10 -3 
at saturation and positive. Lesur & Papaloizou (2010) re- 
port a much lower (of the order of 10 -5 ) and negative an- 
gular momentum transport. This is because of the anelas- 
tic approximation, as the authors themselves point out. 
In that case, the angular momentum " transport" is solely 
due to the 3D instability that taps energy from the vortical 
motion. Compressibility allows for the excitation of spiral 
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Fig. 7. Evolution of enstrophy, kinetic stresses, and vertical ve- 
locities in a 3D baroclinic simulation. The evolution is very sim- 
ilar to the 2D case up to 120 orbits. At that time the vortex goes 
elliptically unstable, and the kinetic energy of vertical motions 
increases by 10 orders of magnitude in less than 10 orbits, but re- 
mains three orders of magnitude lower than the radial rms veloc- 
ity. This 3D elliptical turbulence is very subsonic, and the vortex 
is not destroyed. The level of enstrophy and angular momentum 
transport remain similar to that of a 2D simulation. 
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Fig. 8. Vertical slice of the elliptically unstable vortex core show- 
ing vertical vorticity (left panel) and vertical velocity (right 
panel). The motion in the core constitutes a subsonic turbulence 
at maximum speeds reaching 10% that of sound. 

density waves, which enable positive angular momentum 
transport. 

4.1. Magnetic fields. Interaction of the MRI and the Bl 

The baroclinic instability demonstrated in the past sections 
seems to be able to drive angular momentum transport in 
accretion disks. As such, it could be thought of as an alter- 
native to the MRI. Nevertheless, an important question to 
ask is how the two instabilities interplay. What happens if 
a magnetic field is introduced into the simulation? 

To answer this question, we take a snapshot of 
the quantities at 200 orbits, and add a constant ver- 
tical magnetic field to it, of strength B=5 x 10 -3 
(j8=27~ 1 c^/i^«570). We assume ideal MHD, i.e. perfectly 
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Fig. 9. Angular momentum transport with only the baroclinic in- 
stability in a 3D run and the baroclinic and magneto-rotational 
instabilities, after 200 orbits. The pattern after including the mag- 
netic field is equal to that generated by MRI-only, from which 
we conclude that the BI is irrelevant if magnetic fields are well- 
coupled to the gas. 



coupling of the field to the gas (run P in Table 2). The same 
setup in a barotropic box leads to MRI turbulence with al- 
pha values of the order of a: « 5 x 10 ~ 2 . When the field is 
introduced into the MRI unstable box, the Maxwell stress 
immediately starts to grow, saturating after ^3 orbits, as 
expected from the MRI. The Reynolds stress due to the 
MRI supersedes the stresses due to the BI by one order of 
magnitude (see Fig. 9). The pattern is the same as with an 
MRI-only box. The conclusion is immediate: the BI plays 
little or no role in the angular momentum transport when 
magnetic fields are well-coupled to the gas. This was intu- 
itively expected, since the BI has weak angular momentum 
transport, as well as slow growth rates. The MRI is faster 
by 1 order of magnitude and much stronger. 

In Fig. 10 we show the evolution of energies, enstrophy, 
and temperature, before and after including the magnetic 
field (at 200 orbits). The magnetic energy behaves as ex- 
pected from the MRI, a fast growth and saturation after 
^3 orbits, with most of the energy stored in the azimuthal 
field. The kinetic energy of the turbulence increases by one 
order of magnitude and is more isotropic, also as expected 
from the MRI. The temperature increases by a factor of ~2 
in 15 orbits. This is because the MRI turbulence heated the 
box faster than the thermal relaxation time could bring the 
temperature back to Tq. 

With this experiment we expected to assess the possi- 
bility of synergy between the instabilities, but as far as we 
can tell, none is observed because the MRI alone dictates 
the evolution. 

4.1.1. Vortex destruction by the magnetic field 

If the evolution of the box-averaged quantities brings no 
surprises, the same cannot be said of their spatial distribu- 
tion. In Fig. 11 we plot vorticity at three consecutive orbits 
after insertion of the magnetic field. Magnetic energy is 
also shown. The vortex, which in a nonmagnetic run re- 
tains its coherence indefinitely, is dilacerated when mag- 
netic fields are included. In Fig. 12 we plot ID spatial av- 



erages against time of the vertical enstrophy (upper panel) 
and magnetic energy of the azimuthal field (middle panel). 
A control run where £ = is also shown (lower panel). 
The figure also shows other simulations (discussed later). 
The run in question is the leftmost one, labeled "P". It is 
apparent from the enstrophy plot that the vortex bulges, 
then gets destroyed as the magnetic energy grows. 

To understand this behavior, we examined the state of 
the vortex prior to inserting the field. In Fig. 13 we mea- 
sure the vorticity profile of the vortex. The figure shows 
a slice at the midplane, where we define a box of size SH 
x 2H centered on the vorticity minimum. We used ellipti- 
cal coordinates such that the radius is rv={x x 'c + Vc /x) 1 ^ 2 ' 
where % = a lb is the vortex aspect ratio (a is the semima- 
jor and b the semiminor axis). The coordinates x c and y c 
are rotated by a small angle to account for the off-axis tilt 
of the vortex, (x c ,y c )=R(x — xo,y — yo), where R is the ro- 
tation matrix, and (xo/3/o) are the coordinates of the vortex 
center, found by plotting ellipses and looking for a best fit. 
We find that ^=4 and a rotation of 3° best fits the vortex 
core. Two such ellipses are plotted in the upper panel of 
Fig. 13. 

We then measured the vertical vorticity within the box, 
averaged all vertical measurements for a given radius, and 
box-plotted the z-averaged measurements against ry. The 
box plot uses a bin of Ary = 0.01. The result is seen in 
the lower panel of Fig. 13, with the radius of the ellipses 
drawn in the upper panel. It is seen that the vortex core (in- 
side the inner ellipse) has a vorticity profile that is well ap- 
proximated by a Gaussian, coy = coq exp(— r 2 /2?q), where 
co$ = 0.62, rg = 0.1, and the radii are in elliptical coordi- 
nates. 

We conclude that the vorticity in the core is close to 
uniform (as a Gaussian is very flat near the peak ampli- 
tude). Because the vorticity is finite and close to uniform, 
so is the angular momentum, and thus little radial shear 
should be present in the core. As the MRI feeds on shear, 
one can expect that a patch of constant (or nearly con- 
stant) angular momentum should be stable. Nevertheless, 
examining the vorticity after 2 orbits of the insertion of the 
field, (upper middle panel of Fig. 11), we notice that the 
core did become unstable. This seems to be a signature of 
the magneto-elliptic instability (Lebovitz & Zweibel 2004, 
Mizerski & Bajer 2009), which we consider in the next sec- 
tion. 

4.1.2. Magneto-elliptic instability 

The elliptic instability has been a topic of extensive study 
in fluid mechanics (see review by Kerswell, 2002). First 
studied in the context of absent background rotation 
(Bayly 1986, Pierrehumbert 1986), the effect of the Coriolis 
force was studied by Miyakazi (1993), followed by the 
effect of magnetic fields by Lebovitz & Zweibel (2004). 
The general case, in which both background rotation and 
magnetic fields are present, has recently been studied by 
Mizerski & Bajer (2009). 

These studies have unveiled two regimes of opera- 
tion, which may as well be seen as two different instabili- 
ties. The first de-stabilization mechanism is through reso- 
nances between the frequency of inertial waves and har- 
monics of the vortex turn-over frequency. This instabil- 
ity is three-dimensional, existing for 6 > (the angle 
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Bl+MRI 




Fig. 10. Evolution of box-average quantities (clockwise: kinetic energy, magnetic energy, ens trophy and temperature) before and after 
inserting the magnetic field. The MRI quickly takes over, on its characteristic short timescale. No evidence of synergy between the 
two instabilities is observed. The saturated state of the combined baroclinic+MRI resembles an MRI-only scenario. 




Fig. 11. Evolution of vorticity (upper panels) and magnetic energy (lower panels) in 3D. As the MRI develops, the vortex is destroyed 
by the magnetic field. In a nonmagnetic run, the vortex survives indefinitely. 
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Fig. 12. Time evolution of the ID spatial average of enstrophy (upper panels) and azimuthal magnetic field (middle panels) for the 
runs in Table 2. The lower panels refer to the magnetic field attained in control runs, where £=0. In all runs, the field is seen to grow 
first in the vortex, then in the surrounding flow. This shows that the growth rates of the magneto-elliptic instability are faster than 
those of the MRI. Vortex destruction is apparent in these plots as loss of spatial coherence in the enstrophy plots, and occurred in all 
simulations. The length of the time axis is the same for all simulations, except the control run for run S. 



being the angle between the wavevector of the pertuba- 
tions and that of the vortex motion). Lebovitz & Zweibel 
(2004) show that this instability persists in the presence 
of magnetic fields, and that its effect is twofold. While it 
lowers the growth rates of the elliptically unstable modes, 
the excitation of MHD waves allows for de-stabilization of 
whole new families of resonances. 



The second destabilizing mechanism occurs only when 
the Coriolis force is included (Miyakazi 1993). This in- 
stability is nonresonant in nature and exists only for 0=0 
modes, i.e., oscillations in the same plane of the motion 
of the vortex. Because this plane is associated with a 
"horizontal" (xy) plane (thus k z modes), this destabiliz- 
ing mechanism has been called "horizontal instability". 
As shown by Lesur & Papaloizou (2009), this nonreso- 
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Fig. 13. Vorticity profile of the vortex, prior to the insertion of the magnetic field. We measure the vertical vorticity in the midplane of 
the simulation against the elliptical radius, in the grid points boxed by the thin black line as shown in the upper panel. The modulus 
of the vorticity is plotted in the lower panel. The conclusion is that the vortex core has an angular velocity profile close to uniform, 
with shear where it couples to the Keplerian flow. The dashed lines in the lower panel mark the position of the dotted ellipses in the 
upper panel. They have an aspect ratio ^=4, and mark elliptical radii of ty =0.065 and ry=0.13. The inner one encloses the vortex core. 



nant instability results in exponential drift of epicyclic 
disturbances. It can thus be regarded as an analog of 
the Rayleigh instability, but for elliptical streamlines. For 
a vortex embedded in a Keplerian disk, the modified 
epicylic frequency goes unstable for the range of aspect 
ratios 3/2<^<4. 

Mizerski & Bajer (2009) present the analysis of the gen- 
eral case, including both the Coriolis and Lorentz forces. 
They confirm the previous effects of the Coriolis and 
Lorentz forces in isolation, and find that the horizontal in- 
stability, when present, dominates over all other modes. 
They also find that the magnetic field widens the range of 
existence of the horizontal instability to an unbounded in- 
terval of aspect ratios when 
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is the Rossby number, 

S = (x + X- 1 )/2 
is another measure of the ellipticity, and 



b = 
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(16) 
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(18) 



is a dimensionless parametrization of the magnetic field, 
with k the wavenumber. We can also write b=q /Ro, where 



(19) 



Fig. 14. Shear as the common ground between the magneto- 
rotational and magneto-elliptical instabilities. The distance be- 
tween two points in uniform rotation does not increase if the 
streamlines are circular, i.e., the rotation is rigid (upper figure). 
However, in elliptic streamlines the distance between the two 
points does increase even if the rotation is uniform (lower fig- 
ure). A magnetic field connecting the two points will resist this 
shear, leading to instability depending on the field strength. 



is a more usual dimensionless parametrization of the 
magnetic field, based on the Balbus-Hawley wavelength 
A BH = 2nv A /0 K (Balbus & Hawley 1991, Hawley & 
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Fig. 15. Left. Numerically calculated growth rates of the magneto-elliptic instability for elliptic streamlines of aspect ratio x = 4 as 
a function of the dimensionless magnetic field strength q=kv A /OK and the angle between the wavevector of disturbances and 
the vertical axis. Pure k z modes are the most unstable ones, having a critical wavelength near the predicted ^ cr i t =2y / |Ro|. Weaker 
destabilization exists at intermediate (3D disturbances) for shorter wavelengths. Pure planar disturbances (6=n/2) are stable. 
Right. Growth rates of the k z modes for different aspect ratios. For x = 2 and % = 3 the purely hydrodynamical elliptical instability 
is seen as finite (and high) growth rates as q tends to zero. For x — 4 onwards, the instability is magnetic and has a most unstable 
wavelength near q = 1. The x — 100 curve stands for an approach to the limit of pure shear flow. The growth rate curve calculated 
matches that of the MRI. 



Radial field 



Azimuthal field 



Vertical field 



10" 
10" e 
10" £ 
10" 1C 
10" 12 




200 205 210 215 

t/(2nQr Q ') 



Enstrophy 





- 



200 205 210 215 

t/(2n^) 



Kinetic energy 





200 205 210 215 

t/(2n£lo) 



200 205 210 215 

t/(2na- ') 



200 205 210 215 

t/(2n£lj) 

Run O: Hydro run 
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Strong field, X cr]t >L z , MRI Stable 
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Most unstable wavelength q=\ underresolved 
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Fig. 16. Isolating the vortex magnetic action. The lines show the magnetic runs with the MHD instabilities (magneto-elliptic and 
magneto-rotational) resolved in ideal MHD (P, black solid); unresolved with strong field (Ql, green dot-dot-dot-dashed); most un- 
stable wavelength under-resolved with weak field (R, red dashed); and quenched with resistivity (S, blue dot-dashed). The nonmag- 
netic 3D hydro run is shown as dotted line in the lower panels. Run Ql shows that strong magnetic fields have a stabilizing effect on 
the elliptical instability. Run R shows the magneto-elliptic instability seeming to saturate (at 203 orbits) before the MRI takes over (at 
207 orbits). In runs Q2 and S2, the vortex survives until a channel flow develops in the box. 



Balbus 1991). The analysis of Mizerski & Bajer (2009) also 
assumes that the vortex is of the type 

u x = -Q v (y-y )/x 

u y = O v (x-x )x- (20) 



When there is a magnetic field but the criterion posed 
by Eq. (15) is not fulfilled, Mizerski & Bajer (2009) find that 
the field has an overall stabilizing effect on the resonant 
modes of the classical (hydro) elliptic instability. 
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We can rewrite Eq. (15) in more familiar terms by isolat- 
ing the wavenumber and expressing the criterion in terms 
of Ok and Ro 



0<fc<2|Ro| 1/2 ^. 



(21) 



We estimate the Rossby number of the vortex in Fig. 13. 
Assuming the elliptic flow of Eq. (20), the vorticity is 
coj = 2SOy. The subscript T stands for "total". This dis- 
tinction is necessary because the sheared flow amounts 
to a vorticity of o^box = —3,0^/2. The total vorticity is 
coj = cvv + ^box/ where coy is the vortex's intrinsic vor- 
ticity. By isolating QyS and dividing by Q^, we have 



Ro 



OJy 



(22) 



In the absence of a vortex, the Rossby number is still finite, 
Ro = —3/4, because of the vorticity of the shear flow. In 
this limit, Eq. (21) becomes 

< k < V3Q K /v A , 

which is the criterion for the MRI (Balbus & Hawley 1991). 
As we shall see, the growth rate in this limit also matches 
that of the MRI. This suggests that the MRI is a particular 
case of the magneto-elliptic instability. 

Since we measured coy/Qj^ ~ —0.6, the Rossby num- 
ber is approximately Ro « — 1. These results are compati- 
ble with the Kida solution (Kida 1981) 



n v = 



3Q K 



from which we derive 



LOy 



and thus 
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For #=4, the expressions above yield coy/Q^ = —5/8 = 
—0.625, which matches well the vorticity plateau mea- 
sured in Fig. 13, and Rossby number Ro = —17/16 ~ —1. 
Since the Rossby number is Ro « —1, Eq. (21) implies 
that the horizontal instability in the vortex is present when 
q<2. 

In dimensionless units, we use L z = 0.2, and a resolu- 
tion of N z =128 points in the z-direction, so the wavenum- 
bers present in the box are fcg < k < fc^y, where fcg = 
2rc/ L z = 31 is the largest scale, and fc^y = n/ Az = 2011 is 
the Nyquist scale. The PENCIL CODE needs eight points to 
resolve a wavelength without significant numerical dissi- 
pation, so for practical purposes, the maximum wavenum- 
ber of the inertial range is k^y/4 = 503. Also in dimen- 
sionless units, }io=po=f)K=l, so for Bq=v a =5 x 10 -3 , the 
condition posed by Eq. (21) is k < 400, well within the 
range captured by our box. 

As for growth rates, Mizerski & Bajer (2009) do not un- 
veil a simple expression. The solution has to be computed 



numerically 3 . The growth rate is a function of x> q> an d 
the angle between the z-axis and the wavevector of the 
perturbation. Technically, the Rossby number is also a free 
parameter, but the Kida solution ties the Rossby number 
to the aspect ratio. We present the growth rates f or % = 4 
in the q-6 plane in the left panel of Fig. 15. 

It is seen that the most unstable modes are those of the 
horizontal instability (0 = 0, or k z modes). The right panel 
of Fig. 15 shows the growth rates of these modes for a se- 
ries of aspect ratios. For x = 2 and x — 3, we are in the 
range of existence of the classical (hydro) horizontal insta- 
bility, noted by the fact that fast exponential growth exists 
for v A = 0. For % = 4 onwards, the instability does not 
exist or is too weak in the hydro regime, and the most un- 
stable wavelength is found in the vicinity of q=l. Although 
a critical wavelength exists for k z modes, 3D resonant in- 
stability exists for an unbounded range of wavenumbers, 
albeit with slower growth rates. In the next sections, un- 
less otherwise stated, whenever we mention "magneto- 
elliptic instability" we mean the horizontal, nonresonant, 
magneto-elliptic modes. 

We also calculate the growth rate in the limit of pure 
shear flow, which we approximate numerically by x = 
100. In this case, there are no 3D unstable modes, since 
there is no finite vortex turnover frequency to establish 
resonances. The only instability present is horizontal (fc z ), 
which we also show in the right panel of Fig. 15. The crit- 
ical wavelength is q ~ most unstable wavelength 
q « 1, with growth rate a « 0.75^. One immediately 
recognizes that these properties are the properties of the 
MRI. That the MRI is a limiting case of the magneto-elliptic 
instability makes sense because Eq. (20) with the Kida so- 
lution reduce to a Keplerian sheared flow when % tends 
to infinity. Physically, the destabilization of k z modes of 
the elliptic and magneto-elliptic instability mean expo- 
nential drift of epicyclic disturbances. The equivalent of 
epicyclic disturbances for x —> 00 are radial perturbations 
in a sheared flow. The magneto-rotational and magneto- 
elliptic are in essence the same instability. Another way 
of seeing the common ground between the instabilities is 
by realizing that, although constant angular momentum 
means rigid rotation in circular streamlines, it does not 
mean so when it comes to elliptical streamlines. Figure 14 
illustrates this point. The length of a line connecting two 
points is conserved in uniform circular motion, but not in 
uniform elliptical motion 4 . In other words, uniform ellip- 
tical motion contains shear. A magnetic field connecting 
the two points will resist that shear, leading to instability 
depending on the field strength. 

Judging from the Fig. 15, the growth rates of the 
magneto-elliptic instability at the Balbus-Hawley wave- 
length q = 1 seem to be well-reproduced by a fit 



(7 BH « 0.75O K 



X 



(26) 



3 A script to calculate the growth rates was kindly provided by 
K. Mizerski. 

4 As pointed out by the referee, this is clearly seen when one 
writes the shear stress S s \ x =d x Uy + dyU x and substitutes the Kida 
solution. It yields S sh =— Oy(^ 2 — 1)/X/ which is only zero for 
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i.e., scaled by a factor (x + !)//£ with respect to those of 
the MRI. We hereafter refer to this % = 100 1 curve as 
the MRI limit. 

4.2. Isolating the vortex magnetic action 

As seen in Fig. 15, the wavelength range of the magneto- 
rotational and (horizontal) magneto-elliptic instabilities 
are almost the same for the aspect ratio of interest, leav- 
ing only a narrow range where one instability is captured 
but not the other. However, the growth rates differ, and we 
can explore this fact. The maximum growth rate for % = 4 
is a « 0.95,0^. While the MRI is amplified a millionfold 
in three orbits, the magneto-elliptic instability is amplified 
by more than a billionfold in the same time interval. We 
study in this section limiting cases where the instabilities 
do not grow as fast as in Fig. 11, thus allowing us to bet- 
ter study their behavior. Because the magneto-rotational 
and magneto-elliptic instabilities will both be present in 
the simulations, we loosely refer to them collectively as 
"the MHD instabilities" or just "the instabilities" in the 
next sections. 

4.2.1. Increasing the field strength - Stabilization of elliptic 
instability and channel flows 

We add to the box a vertical field of strength B z =6 x 10 -2 . 
Since the smallest wavenumber of the box is fco=31, we 
have that the critical wavenumber for the MRI is fc/fco = 
0.9 and thus the box is MRI-stable. The critical wavenum- 
ber for the magneto-elliptic instability, according to Eq. 
(21), is fc/fco = 1.13 and thus in principle resolved. We 
aim with this to explore the window between y/3 < q < 
2y|Ro| where the MRI is suppressed but not the horizon- 
tal magneto-elliptic instability. 

We follow the evolution of box-average quantities in 
Fig. 16. The run in question is shown in that plot as green 
dot-dot-dot-dashed lines, and corresponds to run Ql in 
Table 2. 

After insertion of the field, we immediately see a de- 
crease in the box average of the vertical velocities. The ver- 
tical vorticity is unchanged. Radial and azimuthal fields 
decay with the decay of the vertical velocity. A weak 
vertical magnetic field of rms j6=1000 is sustained. Even 
though the analysis provides an elliptical wavelength that 
is shorter than the box length, we do not seem to witness a 
magneto-elliptic instability in operation. In fact, we are in 
the range of stable Rossby numbers, evidence of which is 
that the elliptic instability was stabilized after inserting the 
field. This is not really worrisome considering that the de- 
rived critical wavenumber was so close to fco, an d we made 
some approximations. It is curious, though, that we do 
not see growth in the unstable resonant modes. For wave- 
lengths emcompassing the vortex core (Xy=2H or A x =H/2; 
each of them well-resolved with 32 points), the maximum 
growth rate is at the vicinity of a = 0.330^, yielding a 
millionfold amplification in less than 7 orbits. We ran this 
simulation for 30 orbits after inserting the field. At present, 
we can offer no explanation as to why these modes did not 
become unstable. 

We also test a run with a slightly less strong field, 
B z =3.75 x 10~ 2 (run Q2). In that case, the magneto- 
rotational and magneto-elliptic instabilities have critical 



wavelengths of fc/fcg = 1.47 and fc/fco = 1.78, respectively, 
so both instabilities ought to be resolved. The largest scale 
of the box corresponds to q = 1.18, close to the maximum 
growth rate of both instabilities. The magneto-elliptic in- 
stability has a faster growth rate, so it should be seen first. 

What we witness is quite revealing. The vortex is de- 
stroyed in 4 orbits, while the MRI is still growing in the 
box. A growth of magnetic energy occurs within the vor- 
tex at a very fast pace. The vortex is destroyed when still in 
the phase of linear growth of the instability, owing to the 
development of a conspicuous and strong channel flow 
(Fig. 17). Because the flow in different layers occurs in dif- 
ferent directions, the vortex is stretched apart and loses its 
vertical coherence. 

We notice that, prior to the excitation of the channel 
flow, the elliptical instability in the core was suppressed, 
which is also obvious from comparing the snapshots at 
£=200 and 201 orbits at Fig. 17. 

The run is also shown in Fig. 12 (run Q). We see that 
the growth of magnetic energy occurs earlier in the vortex 
when compared to the surrounding flow, as expected. 

4.2.2. Decreasing the field strength - Vortex MHD turbulence 

Next we checked the behavior of the flow by adding weak 
magnetic fields. The goal was to slow the MHD instabil- 
ities by not resolving their most unstable wavelengths, 
q « 1. Both instabilities thus operate at a slower pace, 
which results in stretching the time interval while one 
(magneto-elliptic in the vortex) is saturating and the other 
(magneto-rotational in the box) still growing. 

The cell size in the z-direction is 1.6 x 10 -3 . We add 
a field of strength B z =1.5 x 10 -3 . The Balbus-Hawley 
wavenumber is /cbh = 667, and thus resolved but within 
the viscous range. The first properly resolved wavenum- 
ber is k « 500, which corresponds to q « 0.75. The run is 
shown in Fig. 16 as dashed red line, and labeled R in Table 
2 and Fig. 12. 

It is seen that the MRI in the Keplerian flow is sup- 
pressed, yet an instability is present. We identify it with the 
magneto-elliptic instability, as it coincides with the vortex 
core going unstable, as shown in the snapshots of Fig. 18. 

The vortex is magneto-elliptic unstable, yet it does not 
seem to lose its spatial coherence. The magnetic field is 
mostly confined to the vortex, which shows as a region of 
high Alfven speeds, when the surrounding Keplerian flow 
is still laminar. The instability is violent, making the vor- 
tex bulge. This is apparent in Fig. 18 as the vortex seems 
to have grown radially from £=203 to £=206 orbits. During 
this period, however, the box average of kinetic energy 
and enstrophy are nearly constant (Fig. 16), so it is not 
clear if this magneto-elliptic turbulence would have led to 
vortex destruction, or if it would have reached a steady 
state. The process just outlined is well-illustrated in Fig. 12 
(run R). One orbit later, the MRI starts to develop in the 
surrounding Keplerian flow (notice the difference between 
these time scales and those of Fig. 11), which corresponds 
to the increase in box-average quantities in Fig. 16 at that 
time. No strong channel flow is excited. The level of vor- 
ticity due to the MRI is nonetheless bigger than that of the 
vortex. The latter eventually becomes inconspicuous in the 
midst of the box turbulence. 
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Fig. 17. The effect of a strong unstable vertical magnetic field in the vorticity column. The field is added at t=200. At first, the effect 
of the field is to stabilize the elliptic turbulence, which is seen in the subsequent snapshots. The disappearance of the vortex at later 
times is caused by the development of a strong channel flow that stretches the column and destroys its vertical coherence. If the 
initial vertical field is stable, the strength of the channel does not grow and the vortex survives indefinitely. 
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Fig. 18. Time series of enstrophy and plasma beta for run R, where the instabilities grow at lower growth rates than in run P (Fig. 11). 
The magneto-elliptic instability grows faster than the MRI, which is seen as the strong turbulence that develops in the core, while the 
underlying Keplerian flow is still laminar. Once the MRI saturates, the strain of its turbulence destroys the vortex spatial coherence. 
It is not conclusive if the vortex would have survived the magneto-elliptic instability had the MRI not destroyed it first. 



We also tested a weaker field, of strength B z =6 x 10 . 
The wavenumbers of the analysis above were then scaled 
by 2.5, so the first resolved wavenumber corresponds to 
q = 0.3. In this case, no significant action was seen. After 
ten orbits, the intensity of the magnetic energy was only 
4 x 10 -9 , accompanied by a merely slight increase in the 
kinetic energy of the vertical velocities (<u^> and <u^> 
remained unchanged). The minimum plasma beta was 
still as high as 10 4 . 



4.2.3. Resistivity 

To test the last case, we used a resistivity high enough that 
the longest unstable wavelength present in the box has a 
magnetic Reynolds number of unity. This wavelength is 
of course L z , the vertical length of the box. The resistivity 
then is such that ReM = L z v A /t] = 1; for a field of strength 
B z =5 x 10 -3 , this magnetic Reynolds number corresponds 
to 7/=10 -3 . This is the same field that was used in the fidu- 
cial MHD run (Fig. 11), of /c B h = 200, so the instabilities 
are resolved in the absence of resistivity. The run is labeled 
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Fig. 19. Time series of enstrophy and plasma beta for runs S, where the instabilities are quenched with resistivity The upper panels 
correspond to a high-resistivity run, where even the longest wavelength of the box is damped. The simulation is similar to a nonmag- 
netized run, the vortex surviving indefinitely In the lower panels we used lower resistivity with Elsasser number A=l. The longest 
wavelength of the box thus has a magnetic Reynolds number of 6, so its growth is not quenched. The magneto-elliptic instability 
grows in the vortex core in a conspicuous k z /ko=2 mode. Part of the field generated is diffused away due to the high resistivity. 
Channel flows eventually develop, destroying the vortex. 



S in Table 2. The results are shown in the upper panels of 
Fig. 19. 

The simulation is not very different from a purely hy- 
dro run. The damped magnetic field only has a slight sta- 
bilizing effect on the elliptical instability. A slight amount 
of the kinetic energy of the core turbulence gets converted 
into magnetic energy, which then diffuses away. The vor- 
tex becomes, at later times, less magnetized than the sur- 
roundings. 

The situation should change when the resistivity is 
lowered slightly, allowing some unstable wavelengths to 
have ReM > 1/ yet still quenching the most unstable 
wavelengths. For that, we set the Elsasser number to A = 
2rcv 2 A / QxT] — 1. The Elsasser number is equivalent to the 
magnetic Reynolds number ReM=^^/ V taking the length 
L as the MRI wavelength, and velocity U as the Alfven ve- 



locity. As such, it is the quantity governing the behavior 
of the MRI (e.g., Pessah 2010). Having A=l corresponds 
to rj=2TCV 2 / Ok/ or V — 1-6 x 10 -4 in dimensionless units. 
The magnetic Reynolds number of the longest wavelength 
is thus « 6. The results are shown in the lower panels of 
Fig. 19. 

The vertical field again has a stabilizing effect on the 
elliptical turbulence. This is seen as a weakening of the 
vertical kinetic energy in Fig. 16, which lasts for two orbits. 
The difference between this run and the more resistive one 
is that thanks to the excitation of magneto-elliptic modes, 
radial and azimuthal fields grow inside the vortex core, 
and a conspicuous k/ko=2 vertical mode appears (lower 
panels of Fig. 19). The field gets looped around the vor- 
tex, initially making the vorticity patch a region of higher 
Alfven speeds. Owing to the high resistivity, however, the 
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Fig. 20. Time series of enstrophy and azimuthal field in the midplane for run S2, of moderate resistivity By action of the magneto- 
elliptic instability, the field initially grows inside the vortex. Due to the resistivity, it then diffuses away from the vortex, coupling to 
the waves excited by it. At later times, exponential growth of the field is seen in the wake. The vortex itself appears unmagnetized. 



field diffuses away (the time for the field to diffuse over a 
scale height is t=H 2 /r/ « 10 orbits). The radial field gets 
sheared into azimuthal by the Keplerian flow. After a few 
orbits, strong magnetic fields are seen in the vortex spi- 
ral waves. At later times, the exponential growth of radial 
and azimuthal fields, as well as the excited z-velocities, are 
seen in these waves. This process is illustrated in Fig. 20. 

A look at the induction equation illustrates the process. 
Under incompressibility and elliptical motion (Eq. (20)), 
the equations for the in-plane field perturbations under 
the influence of a vertical magnetic field and resistivity are 

dB 1 

= B z d z u' x + B' y n v /x + vV 2 B' x (27) 

dB 1 

= B z d z u' y - B' x Q v x + qCt Q B' x + t/V 2 ^. (28) 

The first term in both equations generates field out of 
velocity perturbations. This is the only source term for 
the in-plane field. Under waning velocity perturbations, 
the generation of fields dies out as well. The second term 
is also stretching, but under the vortical motion, thereby 
turning radial fields into azimuthal and vice-versa, at the 
vortex frequency. Its effect is to wrap around the vortex the 
fields generated by the first term. The field then diffuses 
away due to the resistivity. The radial field is sheared into 
azimuthal because of the third term in the azimuthal field 
equation. 

In this simulation, the magnetic Reynolds number of 
the longest wavelength is ReM = ^z^ a /t7=6, so even 
though the most unstable wavelength of the MRI is sup- 
pressed, slower growing wavelengths are present. Since 
they amplify the field, strong channel flows eventually ap- 
pear in the simulation. At 10 orbits the azimuthal field of 
the channel achieves the same strength as that of the field 
in the vortex. The box went turbulent at 15 orbits (Fig. 12), 
slightly after the destruction of the vortex by the channel 
flow. We note, however, that in the control run for the sim- 
ulation in question, the MRI grew slower, only becoming 
noticeable at « 20 orbits (notice the larger range of the 
time axis for the control run). It appears that the field pro- 
duced by the magneto-elliptic instability in the vortex and 
then diffused to the box led to the faster growth compared 
to the £ = control run. 

A simulation where the Reynolds number of the 
longest wavelength was three also shows the same qual- 
itative behavior, albeit in longer timescales. We followed a 
simulation of Reynolds number two for the same time, and 



no growth was seen. The timescale for growth in this case 
may be infinite (stable) or just impractically long. We con- 
clude that resistivity suppresses the magneto-elliptic insta- 
bility when the longest unstable wavelength has a mag- 
netic Reynolds number of the order of unity, as intuitively 
expected. 

4.3. Constant azimuthal field and zero net flux field 

The analysis of the magneto-elliptic instability by Mizerski 
& Bajer (2009) was done for a system thread by a uniform 
constant magnetic field. We seek here to establish the effect 
of a zero-net flux field. As it turns out, the vortex is quite 
unstable to such configurations as well. 

We add a field whose initial value is B=Bq sin(fc z oz)# , 
where fco z =2zr/L z , and Bq = 10 -2 . The run is labeled U 
in Table 2 and Fig. 12. The most unstable wavelength for 
the MRI has k = 100, hence well-resolved. In a barotropic 
box, this field led to saturated turbulence after 3 orbits. The 
critical wavelength for the magneto-elliptic instability has 
k = 200, also well resolved. As shown in Fig. 12, the vortex 
becomes unstable well before the box turbulence starts. In 
1 orbit after insertion of the field, the vortex column has 
already lost its coherence. 

As for an azimuthal field, Ker swell (1994) studied the 
effect of toroidal field on elliptical streamlines, finding 
only a slight stabilizing adjustment of the growth rates of 
the elliptical instability. The analysis, however, only holds 
for the limit of nearly circular streamlines (x —> 1)- Given 
the stark difference in the behavior of vertical fields in dif- 
ferent configurations, there is reason to believe the same 
should apply to azimuthal fields. We add a field B = B§y, 
with Bq = 0.03 (run T in Table 2). Once again, the vortex is 
quickly destroyed, as seen in Fig. 12. 

5. Conclusions 

We model for the first time the evolution of the baroclinic 
instability in 3D including compressibility and magnetic 
fields. We find that the amount of angular momentum 
transport due to the inertial-acoustic waves launched by 
unmagnetized vortices is at the level of oc « 5 x 10 -3 , pos- 
itive, and compatible with the value found in 2D calcula- 
tions. 

When magnetic fields are included and well-coupled 
to the gas, an MHD instability destroys the vortex in short 
timescales. We find that the vortices display a core of 
nearly uniform angular velocity, as claimed in the litera- 
ture (e.g., Klahr & Bodenheimer 2006), so this instability is 
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not the MRI. We identify it with the magneto-elliptic insta- 
bility studied by Lebovitz & Zweibel (2004) and Mizerski 
& Bajer (2009). 

Though Lebovitz & Zweibel (2004) report that the 
magneto-elliptic instability has lower growth rates than 
the MRI, our simulations show the vortex core going un- 
stable faster than the box goes turbulent. That is because 
the presence of background Keplerian rotation allows for 
destabilization of k z modes (horizontal instability), which 
have higher growth rates. We also show that the stability 
criterion and growth rates for the magneto-elliptic insta- 
bility derived by Mizerski & Bajer (2009), when taken in 
the limit of infinite aspect ratio (no vortex) and with shear, 
coincide with those of the MRI. Both instabilities have a 
similar most unstable wavelength, yet the growth rates of 
the magneto-elliptic instability in the range of aspect ratios 
4 < x < 10 are approximately 3 times faster than for the 
MRI. 

After the vortex is destroyed, the saturated state of 
the MRI+BI simulation resembles an MRI-only simulation. 
The same box-averaged values of oc, enstrophy, kinetic, 
and magnetic energies are measured in the two cases. The 
conclusion is that the background entropy gradient plays 
only a small role when magnetic fields are present and 
well-coupled to the gas. The enstrophy produced by the 
BI is four orders of magnitude lower than that produced 
by the MRI. 

We performed a series of numerical experiments to de- 
termine the behavior of the magneto-elliptic instability in 
limiting cases. First, we increased the field so that the crit- 
ical MRI wavelength is bigger than the box. In that case, 
the elliptical turbulence dies out almost immediately after 
inserting the field. We take it as evidence of the stabiliz- 
ing effect of strong magnetic fields on the classical (hydro) 
elliptical instability. When the field is slightly decreased 
so that some unstable wavelengths are resolved, the mag- 
netic field inside the vortex core grows rapidly, leading to 
channel flows that soon break the spatial coherence of the 
vortex column. 

Second, we slow down the instabilities to better study 
the magneto-elliptic instability in isolation. Decreasing the 
growth rate by a factor x stretches the time period between 
their saturations by the same factor. We thus decreased the 
field so that the most unstable wavelength in the box is 
underresolved. In this case, we witness the development 
of magneto-elliptic turbulence in the vortex core only. This 
turbulence was violent, but it is not clear if it would have 
led to destruction of the vortex. After seven orbits, longer 
MRI-unstable wavelengths in the box led to turbulence. 
The vortex was destroyed by the strain of that turbulence, 
bulging and losing coherence, and was eventually lost in 
the turbulent vorticity field of the box. Decreasing the field 
further led to quenching of the magneto-elliptic instability 
as well. 

Third, the instabilities were suppressed with physical 
resistivity, setting the Elsasser number to unity. In this 
case, there is a slight decrease in the kinetic energy of the 
elliptic turbulence, which lasts for two orbits. Meanwhile, 
in-plane magnetic fields develop inside the vortex, loop 
around it, and get diffused away. Vortex destruction hap- 
pens when longer wavelengths in the box, for which 
the magnetic Reynolds number is bigger than one, go 
MRI-unstable. Channel flows develop, and the vortex is 
stretched apart. By increasing the resistivity, the instabili- 



ties are quenched when the longest wavelength of the box 
has a magnetic Reynolds number ReM < 2. 

In addition to uniform vertical fields, we also per- 
formed simulations with azimuthal fields and vertical zero 
net flux fields. These different field configurations also led 
to magneto-elliptic instability in the vortex. 

In view of these results, it is curious that the vortex seen 
in the zero net flux MRI simulations of Fromang & Nelson 
(2005) is stable over hundreds of orbits, a fact left without 
explanation to date. We speculate that this may come from 
a lack of resolution in the global simulation to capture the 
magneto-elliptic modes in the core. 

We conclude that the baroclinic instability is impor- 
tant only when magnetic fields are too weakly coupled 
to the gas. Otherwise they are destroyed by the magneto- 
elliptic instability, channel flows, or by the strain of the 
surrounding MRI turbulence. We thus underscore that our 
results fit neatly in the general picture of the layered ac- 
cretion paradigm in protoplanetary disks. If the MRI su- 
persedes the BI, it thus remains the main source of tur- 
bulence in the active zones where ionization is abundant. 
The active layers are unmodified, whereas the dead zone, 
if baroclinic unstable, is endowed with large-scale vor- 
tices and an associated weak but positive accretion rate of 
a ^ 5 x 10 -3 . This value has to be revised by global sim- 
ulations in view of the limited spatial scale of the shearing 
box. If confirmed, it might be sufficient for a steady state 
to be achieved (Terquem 2008), as long as the borders of 
the dead zone are stable (Lyra et al. 2009). It remains to 
be studied what the conditions are when vertical stratifi- 
cation is included, what the precise transition is between a 
BI dominated dead zone and the MRI-active radii/layers, 
and how the accumulation of solids will proceed inside el- 
lip tically turbulent vortex cores. 
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Appendix A: Linearization of the large-scale 
pressure gradient. 

Given a pressure profile following a power law 



(A.l) 



we linearize it by considering t=Rq + x, and Rq ^> x, as 
below 



V = Po(l-£x/R ). 



(A.2) 



The total pressure includes this time-independent part, 
plus a local fluctuation p' 



p = p + p r 

= p(°) - poZx/R , 

where 



p(°) = p + p', 



(A3) 



(A.4) 



and the superscript "(0)" denotes the part of the pressure 
that has no information on the radial gradient. The pres- 
sure gradient is therefore 



Vp = -po£/Ko* + Vp' = -po^/Rox + Vp(°). (A.5) 



The momentum equation then becomes 



Du 1 
— — = --Vp 



(A.6) 



However, this equation would not work because it is not 
in equilibrium. The second term would continuously inject 
momentum in the box. This is a reflection of the fact that 
the pressure gradient modifies the effective gravity, and 
the rotation curve accordingly. In a global disk, the disk 
settles in sub-Keplerian centrifugal equilibrium. We add 
to the equation a term that embodies this equilibrium 



Du 
~Di 



-Vp 



1 



Po 



Rn 



x. 



(A.7) 



The superscript "(0)" is dropped in Eq. (2). The same 
procedure applies to the energy equation, albeit with a 
small caveat, because the energy equation includes not 
only the pressure gradient but also the pressure itself. In 
adiabatic form, the energy equation is 



dE , 



7 



(7-1) 



p V • u, 



(A.8) 



where E = c v pT is the internal energy. The term on the 
right-hand side is the pdV work. The modifications to the 
energy equation come from this term and from the radial 
advection term. Since E = pi (7 — 1), and p is given by Eq. 
(A.3), we have 

a x E = 3*E(°)-p £/[Ko(7-l)L 

so 



9£(°) 
dt 



+ (m-V)E 



(0) 



(A.9) 



(7 



,(0) 



V • u 



R (7-l) 



We now drop the x-dependent term in the pdV work. 
The main reason is because the term is not periodic, there- 
fore its inclusion in the shearing box can potentially lead 
to long-term effects in the simulation, e.g., uneven heating 
due to the dissipation of waves. It also creates a boundary 
in the box, and waves would be subject to refraction and 
reflection upon reaching it. 

Although the discarding is not physically motivated, 
the loss is not as crucial as it may seem at first. This is be- 
cause the term is not important for the instability cycle; 
that is, the pdV is not acting on the buoyancy itself. Clearly 
a vortical flow in geostrophic balance is divergenceless, 
so the pdV term is zero anyway. Thus, neglecting the x- 
dependent term is not affecting the baroclinic instability. If 
included, the term is only affecting the wave pattern, and 
only slightly. Waves will be subject to more heating (cool- 
ing) upon compression (expansion) on one side of the box 
more than the other. In fact, this could also be expressed 
by a radially varying adiabatic coefficient 7. A test sim- 
ulation was run with the non-periodic term included, in 
order to assess the impact of discarding it. It was seen that 
the term effectively creates, over a long time, a small en- 
tropy jump in the box. Nonetheless, the development of 
the baroclinic instability progressed unhindered, leading 
to no significant different results to vortex formation. Even 
the large-scale sound waves have the same statistical prop- 
erties. Growth rates and saturation levels of the turbulence 
are indistinguishable from those attained in the fiducial 
run. 

We are of course losing consistency for dropping a term 
from the equations. The consistent ways of modeling the 
problem are either by using the Boussinesq approximation 
(as in Lesur & Papaloizou 2010), which filters acoustics out 
and eliminates the pdV term, or by means of global disk 
calculations, where the radial dependencies are only given 
in the initial condition, and are free to evolve in time. This 
is admittedly a limitation of the model, but the benefit of 
the high resolution that can be achieved in the local box, at 
the expense of a term that is not important for the instabil- 
ity cycle, was judged worth the loss in consistency. 

The simplified energy equation reads as 



9£(°) 
dt 



+ {u ■ V) E(°) = -_I_p«>) V • u + 



(7-1) 



R (7-l) 
(A.10) 

i.e., the box has the same temperature all over. The depen- 
dency on u x provides the heating/ cooling that a gas parcel 
would experience in a global model. Because of it, a gas 
parcel heats up when climbing the temperature gradient, 
and cools down when descending it. 

From the definition of entropy (Eq. (7)), we write 



s/c v = In ■ 



and taking the derivative, 



7 In ■ 



Po 



(A.ll) 
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Incompressible shear wave 



1 Ds IDE 



c v Dt 



E Dt 



+ 7V • u 



(A.12) 



where we made use of the continuity equation. 
Multiplying the whole equation by E, we have 



Ds DE 7 _ 

~Dt = ~Dt J^y~^V) U ' 



(A.13) 



so the pdV term cancels when converting the energy equa- 
tion into an equation for entropy. We substitute Eq. (A.10) 
into Eq. (A.13) to find 



Ds(°) 



1 / UxPoZ 



(A.14) 



Dt p T(°) l^o(7 
The superscripts "(0)" are dropped in Eq. (7). 

Appendix B: Testing for aliasing 

One critical task to perform before any attempt to quan- 
tify baroclinic growth of vorticity is to test how well the 
code reproduces the analytical results of shear wave the- 
ory. Compressible and incompressible modes have well- 
defined analytical solutions that may be used to assess the 
presence and quantify the amount of aliasing introduced 
by the scheme. 

Aliasing is a feature of finite-difference codes, which 
occurs when a shear-wave swings from leading (k x < 0) 
to trailing (k x > 0). As the radial wavelength of the wave 
approaches zero, it becomes shorter than the width of a 
grid cell. The signal is lost on the Nyquist scale, and in- 
formation is lost owing to the phase degeneracy that is es- 
tablished. There are an infinite number of possible sinu- 
soids of varying amplitude and phase that are solutions 
(all // aliases ,/ of the correct solution), so spurious power 
can be introduced in the wave. This extra power is then 
transferred from the trailing to the leading wave, which 
will again swing-amplify. 

Because of this, it is possible that aliasing by itself may 
lead to spurious vortex growth. In 2D, the energy spuri- 
ously generated at the aliased swing-amplification has no 
option but to undergo an inverse cascade, the end of which 
is coherent vortices. This is a particular danger for PENCIL, 
because the code is both finite-difference and high-order. 
The high-order nature is in most cases a plus, of course, 
since it leads to little overall numerical dissipation. In this 
case, however, it means that the spuriously added power 
will not be discarded. Lower order codes can be diffusive 
enough that the energy introduced by aliasing may be im- 
mediately dissipated. Indeed, Shen et al. (2006) highlight 
that they do not see any aliasing happening, and suggest 
that this is due to the high degree of numerical diffusiv- 
ity in their code. This is of course a case of two negative 
features canceling each other. 

We combine the best of both worlds by using hyper- 
viscosity. It makes the code dissipative only where it is 
needed (on the grid scale), with the extra benefit of be- 
ing able to control how much dissipation is added to the 
solution. We note that Oishi & Mac Low (2009), also us- 
ing the Pencil code, tested the numerical evolution of 
incompressible 2D as well as compressive 2D and 3D dis- 
turbances against their analytical solutions. They found 
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Fig. B.l. Testing the numerical scheme for aliasing, a feature of 
finite difference methods, which spuriously increases the energy 
of a wave that swings from leading to trailing. At the resolutions 
and mesh-Reynolds numbers we use, this effect is successfully 
suppressed. The wave has a wavenumber k X/ o=— 16n/L x , which 
means 4, 8, 16, and 32 points per wavelength at resolutions 32, 
64, 128, and 256, respectively. The wiggling is due to excitation of 
compressible modes not present in the analytical solution. 



aliasing unimportant when using hyperviscosity. We re- 
peat here the incompressible test (showed in Fig. B.l), and 
refer to Oishi & Mac Low (2009) and Johansen et al. (2009) 
for the other tests. 

The analytical solution of the incompressible shear 
wave is (Johnson & Gammie 2005) 



Su x Su Xf0 ^ 



(B.l) 



where k 2 =k 2 x + k 2 , fcy=const, and 

k x = k x (t) =k X/0 + qQ kyt. 
The condition of incompressibility kjSuj 



(B.2) 
thus 

dictates that Suy=—k x Su x /ky. The solenoidality of 
the wave is guaranteed by initializing the velocity 
field through a streamf unction w=V x (ipz), with 
ip=Ak~ r sin (k x x + k y y). We use the setup of Shen et 
al. (2006), k Xr0 =-16n/E x , ky=4n/Ey, with L x =L y =0.5, 
A=10" 4 , c s =n=p =l, y=7/5, and ^=3/2. 

We follow the evolution of the x-velocity in a point of 
the grid, and plot the result in Fig. B.l, checking for differ- 
ences due to resolution and initial mesh Reynolds number 
Re, defined by the hyperviscosity coefficient 



Re 



Su X/ q dx 5 

1/(3) 



(B.3) 



We see that no aliasing is detected at Re=0.01. The sig- 
nal is viscously damped before it can be swing-amplified. 
At Re=0.1, aliasing occurs at resolution 32 2 . However, the 
added signal suffers a severe (hyper-)viscous damping, so 
that the next swing amplification has little left to work 
with. When increasing the Reynolds number to Re=l, 
aliasing now happens also for resolution 64 2 , but the solu- 
tion is decaying. Aliasing is not detected at any Reynolds 
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number for resolution 256 2 . A run with Re=oo for resolu- 
tion 32 2 (not showed) kept periodic intervals of aliasing 
without dissipation but without net growth, at least until 
£=1000. 

The panels in Fig. B.l are for £ =0. When we use non- 
zero £, the aliased solutions show an increase in ampli- 
tude. The same does not happen for the solutions at low 
Reynolds numbers. The wiggling seen at higher resolution 
in Fig. B.l is accompanied by changes in the density, which 
leads us to conclude that they come from the excitation of 
compressible modes not present in the analytical solution. 
We use a resolution of 256 2 for 2D runs. At this resolution, 
we rest assured that aliasing is not happening for the con- 
sidered Reynolds numbers. 



